#include <iostream>
#include <cmath>

#include "CScipyDistance.h"

const float eps = 1e-5;

void ScipyDistance::countNotBool(VectorXf u, VectorXf v, float &tt, float &tf, float &ff, float &ft){
    ArrayXf u_arr = u.array();
    ArrayXf v_arr = v.array();
    ArrayXf not_u = 1.0f - u_arr;
    ArrayXf not_v = 1.0f - v_arr;

    ff = (not_u * not_v).sum();
    tf = (u_arr * not_v).sum();
    tt = (u_arr * v_arr).sum();
    ft = (not_u * v_arr).sum();
}

float ScipyDistance::yule(VectorXf u,  VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: yule distance input should have same nonzero size\n");
        return 0;
    }

    float tt, ff, tf, ft;
    countNotBool(u, v, tt, tf, ff, ft);
    return 2.0f * tf * ft / (tt * ff + tf * ft + eps);
}

float ScipyDistance::braycurtis(VectorXf u, VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: braycurtis distance input should have same nonzero size\n");
        return 0;
    }

    ArrayXf u_arr = u.array();
    ArrayXf v_arr = v.array();

    return (u_arr - v_arr).abs().sum() / ((u_arr + v_arr).abs().sum() + eps);
}

float ScipyDistance::sokalmichener(VectorXf u, VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: sokalmichener distance input should have same nonzero size\n");
        return 0;
    }

    float tt, ff, tf, ft;
    countNotBool(u, v, tt, tf, ff, ft);
    return 2.0f * (tf + ft) / (tt + ff + 2.0f * (tf + ft) + eps);
}

float ScipyDistance::kulsinski(VectorXf u, VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: kulsinski distance input should have same nonzero size\n");
        return 0;
    }
    float tt, ff, tf, ft;
    countNotBool(u, v, tt, tf, ff, ft);
    int n = u.size();
    return (tf + ft - tt + n) / (tf + ft + n + eps);
}

float ScipyDistance::russellrao(VectorXf u, VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: russellrao distance input should have same nonzero size\n");
        return 0;
    }

    ArrayXf u_arr = u.array();
    ArrayXf v_arr = v.array();
    float tt = (u_arr * v_arr).sum();
    float n = float(u.size());
    //printf(" vars: %f, %f", tt, n);
    return (n - tt) / (n + eps);
}

float ScipyDistance::dice(VectorXf u, VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: dice distance input should have same nonzero size\n");
        return 0;
    }

    float tt, ff, tf, ft;
    countNotBool(u, v, tt, tf, ff, ft);
    return (tf + ft) / (2.0f * tt + tf + ft + eps);
}

float ScipyDistance::canberra(VectorXf u, VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: canberra distance input should have same nonzero size\n");
        return 0;
    }

    ArrayXf u_arr = u.array();
    ArrayXf v_arr = v.array();

    return ((u_arr - v_arr).abs() / (u_arr.abs() + v_arr.abs() + eps)).sum();
}

float ScipyDistance::chebyshev(VectorXf u, VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: chebyshev distance input should have same nonzero size\n");
        return 0;
    }

    return (u.array() - v.array()).abs().maxCoeff();
}



float ScipyDistance::moment(const VectorXf &u, int n){
    if(n==0) return 1.0f;
    if(n==1) return 0.0f;
    return (u.array() - u.mean()).pow(n).mean();
}

float ScipyDistance::skew(VectorXf u){
    float m2 = moment(u, 2);
    float m3 = moment(u, 3);
    if(abs(m2) < eps) return 0.0f;
    return m3 / pow(m2, 1.5);
}

float ScipyDistance::kurtosis(VectorXf u){
    float m2 = moment(u, 2);
    float m4 = moment(u, 4);
    if(abs(m2)<eps) return -3.0f;
    return m4 / pow(m2, 2) - 3.0f;
}

float ScipyDistance::jaccard(VectorXf u, VectorXf v){
    if(u.size() != v.size() || u.size() == 0 || v.size() == 0){
        printf("ERROR: jaccard distance input should have same nonzero size\n");
        return 0;
    }

    float tt, ff, tf, ft;
    countNotBool(u, v, tt, tf, ff, ft);
    return 1.0f * (tf + ft) / (tt + ft + tf + eps);
}
